OUTP-97-60-P 



|hep-ph/9804285 



Extremely high energy cosmic rays from relic particle decays 



Michael Birkel & Subir Sarkar* 

Theoretical Physics, University of Oxford, 
1 Keble Road, Oxford 0X1 3NP, UK 



Abstract 



The expected proton and neutrino fluxes from decays of massive metastable 
relic particles are calculated using the HERWIG QCD event generator. The 
predicted proton spectrum can account for the observed flux of extremely high 
energy cosmic rays beyond the Greisen-Zatsepin-Kuzmin cutoff, for a decay- 
ing particle mass of C(10 12 ) GeV. The lifetime required is of O(10 20 ) yr if such 
particles constitute all of the dark matter (with a proportionally shorter life- 
time for a smaller contribution). Such values are plausible if the metastable 
particles are hadron-like bound states from the hidden sector of supersym- 
metry breaking which decay through non-renormalizable interactions. The 
expected ratio of the proton to neutrino flux is given as a diagonistic of the 
decaying particle model for the forthcoming Pierre Auger Project. 
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I. INTRODUCTION 



It has been known for some time that interactions on the 2.73 K blackbody cosmic microwave 
background (CMB) will severely degrade the energies of cosmic ray nucleons with energies beyond 
~ 5 x 10 19 eV — the Greisen-Zatsepin-Kuzmin (GZK) cutoff [|J. It was therefore very surprising 
when the Fly's Eye atmospheric fluorescence detector reported the observation of an extremely 
high energy cosmic ray (EHECR) event with an energy of (3.0 ± 0.9) x 10 20 eV 0. This was 
followed by the detection of a (1.7 - 2.6) x 10 20 eV event by the AGASA air shower array H . These 
discoveries substantiated earlier claims from the Volcano Ranch Q, Haverah Park || and Yakutsk 
H air shower arrays that cosmic rays do exist beyond the GZK cutoff. About a dozen such events 
are now known. Detailed accounts of the data may be found in recent reviews [0]. 

In Figure [l] we show the EHECR spectrum for energies exceeding 10 18 eV note that the fluxes 
have been multiplied by E s . It is believed that cosmic rays with energies up to ~ 5 x 10 18 eV, the 
so-called 'ankle', are predominantly of galactic origin, possibly accelerated by the Fermi mechanism 
in supernova remnants ||. Above this energy, the spectrum flattens and the composition changes 
from being mostly heavy nuclei to mostly protons. Such a correlated change in the spectrum 
and composition was first established by the Fly's Eye experiment [|| and Figure |l| shows their 
suggested two-component fit to the data. The new component which dominates at energies beyond 
~ 5 x 10 18 eV is isotropic and therefore cannot possibly originate in the galactic disk |IT||TT| ]. 
However it also extends well beyond the GZK cutoff raising serious problems for hypothetical 
extragalactic sources. Because of the rapid energy degradation at these energies through photo- 
pion production on the CMB, such sources must exist within ~ 500 Mpc, in fact within ~ 50Mpc 
for the highest energy Fly's Eye event [12|. For heavy nuclei, the energy loss is less severe according 
to a revised calculation so the range may extend upto ~ 100 Mpc. General arguments p^Jl5| 
provide correlated constraints on the magnetic field strength and spatial extent of the region 
necessary to accelerate particles to such high energies and these requirements are barely met by 
likely astrophysical sites such as active galactic nuclei and the 'hot spots' of radio galaxies fi~6| . 
Moreover there are few such sources close to us and no definite correlations have been found 
between their locations and the arrival directions of the most energetic events fl7|]Io| . It has been 
speculated that gamma-ray bursts which too are isotropically distributed, may be responsible for 
EHECRs |l8|]. However since these are at cosmological distances, one would expect to see the GZK 
cutoff in the cosmic ray spectrum contrary to observations (cf. ref. [|19|j). 

Some of the above arguments may be evaded if the EHECR events are due not to nucleons 
but neutral particles such as photons and neutrinos. Although high energy photons also suffer 
energy losses in traversing the CMB and the extragalactic radio background, there is no threshold 
effect which would cause a cutoff near the GZK value |2Q| . However the observed shower profile 
of the highest energy Fly's Eye event Q argues against the primary being a photon since it 
would have interacted on the geomagnetic field and started cascading well before entering the 
atmosphere [21|. The observed events are also unlikely to be initiated by neutrinos as they all have 
incident angles of less than 40° from the zenith and thus too small a path length in the atmosphere 
for interactions |J22[|. This argument may be evaded if neutrinos become strongly interacting at 
high energies due to new physics beyond the Standard Model [ 23 , 24 1 , but such proposals are 
found not to be phenomenologically viable |25| (although this is disputed f2(|). (Alternatively, 
the propagating high energy neutrinos could annihilate on the relic cosmic neutrino background, 
assumed to have a small mass of O(0. 1) eV, to make hadronic jets within the GZK zone f27|.) 
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Other exotic possibilities have been suggested, e.g. monopoles [28], stable supersymmetric hadrons 
[29] and loops of superconducting cosmic string ('vortons') |30|. However these possibilities have 
many phenomenological problems [31,32] and we do not discuss them further. 

Thus one is encouraged to seek 'top-down' explanations for EHECRs in which they originate 
from the decay of massive particles, rather than being accelerated up from low energies. The most 
discussed models in this connection are based on the annihilation or collapse of topological defects 
such as cosmic strings or monopoles formed in the early universe [33|-36]. When topological defects 
are destroyed their energy is released as massive gauge and Higgs bosons which are expected to have 
masses of O(10 16 ) GeV if such defects have formed at a GUT-symmetry breaking phase transition. 
The decays of such particles can generate cascades of high energy nucleons, 7-rays and neutrinos. 
A more recent suggestion is that EHECRs arise from the decays of metastable particles with masses 
mx ~ 10 13 — 10 16 GeV which constitute a fraction of the dark matter p7| . These authors suggest 
that such particles can be produced during reheating following inflation or through the decay of 
hybrid topological defects such as monopoles connected by strings, or walls bounded by strings. 
The required metastability of the particle is ensured by an unspecified discrete symmetry which is 
violated by quantum gravity (wormhole) effects. Another suggestion is that the long lifetime is due 
to non-perturbative instanton effects |38]]. In ref. a candidate metastable particle is identified 
in a 5*7(15) GUT. 

A generic feature of these 'top-down' models is that the EHECR spectrum resulting from the 
decay cascade is essentially determined by particle physics considerations. Of course the subsequent 
propagation effects have astrophysical uncertainties but since the decays must occur relatively 
locally in order to evade the GZK cutoff [37], they are relatively unimportant. Thus although the 
proposal is speculative, it is possible, in principle, to make reliable calculations to confront with 
data. In this work we consider another possible candidate for a relic metastable massive particle 
[40] whose decays can give rise to the observed highest energy cosmic rays. First we discuss (§ ||) 
why this candidate, which arises from the hidden sector of supersymmetry breaking, is perhaps 
physically better motivated than the other suggested relics. We then undertake (§ |T|) a detailed 
calculation of the decay cascade using a Monte Carlo event generator to simulate non-perturbative 
QCD effects. This allows us to obtain a more reliable estimate of the cosmic ray spectrum than has 
been possible in earlier work on both topological defect models [34] and a decaying particle model 
J37j. We confront our results with observations and identify the mass and abundance/lifetime 
required to fit the data. We conclude (§ with a summary of experimental tests of the decaying 
particle hypothesis. 



II. MASSIVE, METASTABLE DARK MATTER FROM THE HIDDEN SECTOR 

Soon after the discovery of the anomaly-free heterotic superstring theory in ten dimensions 
based on the gauge group Eg x E$, it was pointed out Q] that in the physical low energy theory 
(where a grand unified Eq or O(10) group is broken by Wilson lines), the minimum value of 
magnetic charge is not the Dirac quantum 2-7r/e but an integral multiple thereof. Conversely, the 
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minimum electric charge is smaller than the electron charge e by the same ratio. [] This was found 



to be a generic feature of all superstring models based on a level-one Kac-Moody algebra [42]. In 



view of the severe experimental upper bounds on the relic abundance of fractional charges [43], 



this posed a potential embarrassment for superstring phenomenology [44]. 

A simple solution to the problem of fractional charges (with an obvious historical analogue 
in quarks and QCD) is to confine them and it was shown that this can be done in the hidden 



sector of supersymmetry breaking in the framework of the SU(5) (g> U(l) unification model [45]. 
In this model, all fractionally charged states have charges |Q e m| = \ and are placed in 4 or 6 
representations of a hidden SU(4:) gauge group which becomes strong at a scale A4 ~ 10 12 GeV and 
in 10 representations of a hidden 50(10) group which becomes strong at a scale A10 ~ 10 15 GeV. 
This results in integer-charged particles — 'cryptons' — which may be 2-constituent mesons, 3- 
constituent baryons or 4-constituent 'tetrons' [|46|]. Some of these mesons could be light (in analogy 
to the pion of QCD) but most of the states should be heavy with masses of order the confinement 
scale A. (Other possibilities for stable superstring relics have been discussed in ref. pTfl .) 

The constituent fields have very few renormalizable (N = 3) superpotential interactions, so 
most of these states can only decay via higher-order (N > 4) superpotential terms. Generically, 
crypton lifetimes are expected to be 

1 /M\ 2 ^" 3 ) 

t x ^ , mx~A, (1) 



m x \m x , 

where, M = Mp/\/87r ~ 2.4 x 10 18 GeV is the normalized Planck scale, giving 

t 4 ~ l ( 12 ^- 8 °) yr, r 10 ~ 10( 6JV "'- 65 ) yr, (2) 

for SU{A) and SO(W) bound states respectively. Thus r 4 > 1 sec(10 16 yr) for N 4 > 6(8) and 
no ^ 1 sec(10 16 yr) for A^o > 10(14). Detailed studies of the possible effects of decays of relic 



cryptons on primordial nucleosynthesis and the CMB spectrum [48], as well as on the diffuse 7-ray 



background [48|,49| have established that such particles, if they survive as relics of the Big Bang, 
must either decay well before nucleosynthesis or have lifetimes longer than the age of the universe 
(t ~ 10 10 yr). In the latter case, if such particles make an interesting contribution to the dark 



matter, their lifetime is further required to exceed ~ 10 16 yr in order to respect experimental bounds 



on the flux of high energy neutrinos expected from their decays [48,50]. It is seen from eq.(g) that 
these constraints favour SU(4) mesons over their 50(10) counterparts as possible constituents of 
the dark matter. It is then natural to contemplate the possibility that such cryptons with a mass 
of 771,4 ~ 10 12 GeV and a lifetime r 4 > 10 16 yr are also responsible for the observed highest energy 
cosmic rays. 

Recently the above discussion has been extended to other massive metastable particle can- 



didates in superstring/M-theory [40|. These authors discuss constructions with higher-level Kac- 



Moody algebras (necessary to accommodate adjoint Higgs representations in (unified) models other 



x The vacuum state of a physical theory in this scheme must be M 4 x K where M 4 is four- 
dimensional Minkowski space and K is some compactified six-manifold. Such fractionally charged 
states exist because K is not simply connected — these are states in which a closed string wraps 
around a non-contractible loop in K. 
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than SU(5) <8> ^(1)) and note that similar metastable bound states occur in such models. They go 
on to consider other candidate particles in M-theory such as Kaluza-Klein states associated with 
extra dimensions but find that these are not as attractive, being either too heavy or too unstable. 
They suggest that although the SU(5) x U(l) model [^5| discussed above was constructed in the 
weak coupling limit, it may be elevated to an M-theory model in the strong coupling limit. The 
SU(4) tetrons are then still the most likely candidates for massive metastable dark matter with 
the modification that the Planck scale M in Eq.(|l|) may be replaced by a somewhat smaller scale. 

The main reason why this possibility was not seriously entertained earlier concerns the ex- 
pected relic abundance of such massive particles. If cryptons were maintained in chemical equi- 
librium in the early universe through self-annihilations, their present energy density is given by 
the usual 'freeze-out' calculation as inversely proportional to the (velocity-averaged) annihilation 
cross-section [51]. Estimating this to be (a ann v) ~ m^ 2 we see that equilibrium would have been 



established if the annihilation rate exceeded the Hubble expansion rate (H ~ T 2 /M), i.e. at 
temperatures 

T > Tdcc ~ HM/mxY (3) 
The relic abundance is then simply estimated as the equilibrium value at decoupling: 

mx y ( 4 ) 

10 12 GeVy W 
This is the basis for the conclusion that no stable relic particle may have a mass in excess of 



~ 10 5 GeV without 'overdosing' the universe, i.e. contributing fix > 1 ||51,47|. This does not 
necessarily apply to cryptons since a period of inflation should have diluted their abundance to 
essentially zero, along with monopoles and other such supermassive relics. If the reheating tem- 
perature following inflation is restricted to be Tr < 10 9 — 10 10 GeV in order not to produce too 



many gravitinos [52,53|, cryptons would not have been generated afterwards. 



However it has been recently recognized that in supersymmetric cosmology, there is likely to be 
a late stage of 'thermal inflation' [54] due to symmetry breaking along flat directions at intermediate 
scales |55|,|56]]. |^| This would adequately dilute the abundance of thermally generated gravitinos 
following inflation so the bound quoted above on the reheating temperature is no longer valid and 
the value of Tr may be much higher. In that case cryptons even as massive as 10 12 GeV may 



2 This was initially considered to be an 'entropy crisis' |3^] since it would dilute any baryon 
asymmetry generated at the GUT scale. However there are now several plausible mechanisms for 



low temperature baryogenesis [57,58] which may operate after thermal inflation. 



3 The vacuum energy V{4>) of the scalar 'inflaton' field is constrained to be F 1 / 4 /e 1 / 4 ~ 2.7 x 1CT 2 M 
by the anisotropy in the CMB observed by COBE, where the slope parameter e = (M 2 /2)(V' /V) 2 
is required to be <C 1 to permit inflation to occur [ |59|] . (The number of e- folds of expansion until 
the end of inflation is just N = J^ cnd d(f>/M\/2e and this should exceed ~ 50 — 60 in order to solve 
the flatness and homogeneity problems of the standard cosmology.) The reheat temperature Tr 
can, in principle, have been as high as ~ y 1 / 4 although it is usually considerably smaller since the 
inflaton field is very weakly coupled in most inflationary models. 
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well have been brought back into thermal equilibrium during reheating after inflation and survived 
with the huge relic abundance (||). However thermal inflation would also have diluted this to an 
acceptable level as was noted in ref. [ SO | ; to obtain fix ~ 1) the number of e- folds of thermal 



inflation required is just 

iV therm = ~ln(10 14 )~ll. (5) 

This fits in well with the expectation that A^herm ~ ^ln(E/mw) ~ 10 — 15 for the intermediate 
scale S in the range (10 -7 — 10 _2 )M |54j. Of course given the uncertainty in the value of A^herm 
(and indeed the possibility that there may be more than one such epoch), fix could well have been 
reduced to a negligibly small value. 

Another possibility is that massive particles such as cryptons were never in thermal equilibrium 
but were created with a cosmologically interesting abundance due to the varying gravitational field 
during (primordial) inflation [^J6r|]. A cosmologically interesting relic abundance then arises for 
mx ~ (0.04 — 2)H where H ~ 10 11 — 10 14 GeV is the likely Hubble parameter during inflation 
plf . This is certainly very encouraging but it should be remembered that a later stage of thermal 
inflation would dilute such an abundance to a negligible level, as discussed above. 

It is clear that the relic abundance of massive particles such as cryptons will necessarily be 
very uncertain given our ignorance of the thermal history of the universe prior to nucleosynthe- 
sis. However as the above discussion illustrates, there are two complementary ways in which a 
cosmologically interesting abundance may result so we may reasonably consider such particles as 
candidates for the dark matter. We now move on to discuss whether relic cryptons can indeed 
be the source of the EHECR by determining the expected spectrum of high energy particles from 
their decays. 



III. COSMIC RAYS FROM MASSIVE PARTICLE DECAY 

To calculate the expected flux of cosmic rays from decays of massive particles such as cryptons, 
we must consider the contribution from both decaying particles in the halo of our Galaxy as well as 
those elsewhere in the universe. Since such massive particles would behave as cold dark matter and 
cluster efficiently in all gravitational potential wells, their abundance in our galactic halo would be 
enhanced above their cosmological abundance by a factor 

/ cos ee ntVnT ■ (6) 

Note that fix = mx n x s / Pent where p cn t — 1.054 x 10 _4 /i 2 GeV cm~ 3 is the critical density in 
terms of the present Hubble parameter h = //o/lOOkmsec -1 Mpc~ 4 . If for simplicity we assume 
a spherical halo of uniform density, 

-Rhalo ~ lOOkpc, p hal ° ~ 0.3 GeV cm" 3 , (7) 

then / cos ~ 3 x 10 4 /i -2 and the number density of cryptons in the halo is 

^x^cW^W^r. (B) 
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The actual density of dark matter in the halo must of course fall off as r~ 2 to account for the flat 
rotation curve of the disk but we do not consider it necessary at this stage to investigate realistic 
mass models. Thus the universal density of cryptons is smaller than the halo density by about 
the same numerical factor by which the distance to the horizon (~ 3000/i _1 Mpc) exceeds the halo 
radius, so the extragalactic contribution to the EHECR flux from decaying cryptons cannot exceed 



the halo contribution. In particular, the GZK cutoff scale for protons [12| or heavy nuclei [13] are 
all much smaller than the horizon distance, so only the halo contribution need be considered, as 
was emphasized in ref. |57|]. Only for neutrinos would the extragalactic component be comparable 



in magnitude [5C[. Henceforth we restrict ourselves to considering crypton decays in the halo alone. 
Now the injection spectrum from particle decay is, to a good approximation, 



dn x 



dt 



m 

dE 



n 



halo 
X 



2 dNj 



for lifetimes longer than the age of the universe (tx S> to). Here 

E 2E 



jet 



mx 



(9) 



(10) 



is a measure of particle energy (assuming 2-body decays) and the fragmentation function diVj/dx 
is the average number of particles i released per decay, per unit interval of x, at the value x. The 
flux at Earth is then 



ji{E) 



47T 



Rh^i(E). 



(11) 



The final state particles which interest us most are 'protons' and neutrinos/antineutrinos where 
the former includes other nucleons, e.g. antiprotons and neutrons, since they all interact similarly 
in the Earth's atmosphere. To compare with observations we multiply the fluxes by E s and define 



I p (E)=j p (E)E 3 -- 
I V (E) = j u (E)E 3 -- 



1 n 



halo 
X 



4tt tx 



-R 



halo 



1 n 



halo 
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4tt tx 



-R\, 



alo 



2 dN p 
mx dx 

2 dN v 
mx dx 



-E 6 



E' A 



(12) 



For photons and electrons/positrons, propagation energy losses are substantial even within the halo 
and we do not attempt to determine these. However their injection spectra from particle decay are 
given by the same computation as for protons and neutrinos, to which we now turn. 



A. Computing the fragmentation functions 

Heavy particles, whether GUT-scale bosons (in topological defect models), cryptons or other 
hypothetical massive particles, will decay into quarks and leptons. The quarks will hadronize 
producing jets of mostly pions with a small admixture of nucleons and antinucleons. The neutral 
pions will decay to give photons while charged pion decays will yield neutrinos and antineutrinos 
in addition to leptons. Thus the final spectrum of the decay produced particles will be essentially 
determined by the 'fragmentation' of quarks/gluons into hadrons. This is a non-perturbative QCD 
process and it has not been possible to calculate it by analytic means. Usually phenomenologically 
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motivated approximations are used to model experimental data on inclusive jet multiplicities and 
scaling violations |63| . 

So far, authors of proposals involving heavy particle decay, e.g. in the context of topological 
defect models [ 



have employed a hadronic fragmentation function suggested by Hill [ 33 1 



dN, 



(Hill) 



dx 



0.08 



exp[2.6Vln(l/x)](l 



(13) 



It is further assumed that 3% of the hadronic jets from massive relic particle decays turn into 
nucleons, while the other 97% are pions which decay into photons and neutrinos. This was based 



on the leading logarithm approximation of QCD [64] applied to experimental data from PETRA 



on jet production in e + e" collisions at tens of GeV. The estimated jet multiplicity from gluon 
fragmentation was convoluted with the gluon distribution to determine the total hadron yield; to 
estimate the spectrum, it was assumed that the first moment of the distribution is normalized 
to unity and the large x behaviour was guessed to be (1 - x) 2 @. As we shall see, the Hill 
fragmentation function (|i~3| ) significantly overestimates the yield of high x final states from the 
decay of very massive particles and, moreover, photons and neutrinos are actually produced with 
a spectrum quite different from that of nucleons. Thus the decay spectra derived using eq.([l^) for 
topological defect models [ 34 1 are inaccurate. 

Subsequently, another form called the Modified Leading Logarithm Approximation (MLLA) 
which gives a better description of data at low x has been proposed |33) ; a gaussian approximation 
to this is 



dN, 



(MLLA) 



h 



K N 



where Kn is a constant and 



2a z 



dx 



7T 



• cxp 



ln 2 (x/x E 
2^" 



(14) 



2tt 



21 V 3o§(«) 



0.09 



In 



m 



X 



A 2 



3/2 



(15) 



with x max = \/ A/mx and A = 0.234 GeV. This fragmentation function is employed by the authors 
of ref. [37 1 to compute the spectrum from relic particle decays; they determine Kn by requiring 



that the integral of x d jV^ MLLA ^ /d x over the range x G [0, 1] be equal to the fraction of the energy 
transferred to hadrons. However this procedure is not exact as the form ( |i~4|) is inapplicable for 
large x and therefore cannot be normalized in this manner. Thus the shape of the cosmic ray 
spectrum computed |J7[] by this method for decaying particles is only reliable for small x and its 
normalization uncertain. 

Given the importance of determining the energy spectrum accurately, we decided to improve 
on these approximate formulations by using the standard tool employed by experimental high 
energy physicists to study QCD fragmentation, viz. a Monte Carlo event generator. Here the non- 
perturbative hadronization process is simulated on a computer by a well tested phenomenological 
model |66|. Although this requires extensive numerical calculations, it is the only means by which 
successful contact can be made between theory and experimental data. We chose the programme 
HERWIG pBl (Hadron Emission Reactions With Interfering Gluons) which incorporates the cluster 



model for hadronization and is based on a shower algorithm of the 'jet calculus' type 66]. To check 



8 



our results we also ran the JETSET programme [S7| and found good agreement over the energy 
range where comparison was possible. However for the very high energies studied in this work, 
HERWIG proved to be more suitable for reasons of computing time [68]. Even so the calculations 
described here took several months on a Digital Alpha workstation. 

For definiteness, we assume the heavy particles to decay into a quark-antiquark pair with unit 
branching ratio. The quark and antiquark, each carrying away energy mx/2, form jets which lead 
to the generation of many particles through cascading, hadronization and decays of some of the 
generated particles. This can be simulated by HERWIG via the annihilation process e + e~ — * qq 
with center-of-mass energy = mj, where q stands for all six kinematically allowed quark 
flavours. The event generator outputs kinematical details of all final state particles, e.g. protons, 
photons and leptons (electrons, positrons and neutrinos). We divided the x-range into 100 bins of 
width Ax = 0.01. After each event simulation the number of protons, neutrinos, photons as well 
as electrons and positrons per energy bin was counted. We ran 10000 events for each of the masses 
mx = 10 3 , 10 5 , 10 7 , 10 9 and 10 11 GeV. After all events had been run, the particle numbers in the 
bins are divided by the bin width and the number of events, in order to obtain the fragmentation 
functions d Ni/dx. Apart from altering some relevant parameters in the computer code to allow it 
to run at the high energies studied here, we also switched off initial state radiation since it is not 
relevant for the present study. Unfortunately, it was not feasible to study the high x behaviour of 
the fragmentation functions for decaying particle masses higher than 10 11 GeV because of numerical 
convergence problems in the computer code. (Already for masses exceeding 10 9 GeV quadruple 
precision had to be used.) Hence we have had to extrapolate the fragmentation functions to high 
x for very heavy masses as described later. 

First we show the proton fragmentation function obtained from the HERWIG runs in Figure ^ 
to illustrate that it depends on the decaying particle mass, contrary to the approximation (|l3| ) 
employed in previous work on topological defects p4|]. Rather than being constant, it decreases 
with increasing mx for x > 0.1, while at very low x it increases with increasing mx- The large 
fluctuations at x > 0.5 are due to the fact that relatively few particles are produced with such high 
energies despite the 10000 events per simulation. We note also that the shape differs significantly 
at high x from the approximation used in ref. |^J. 

In Figure |3] the fragmentation functions for protons, photons, neutrinos and electrons are 
compared for mx = 10 11 GeV. It is seen that at very low x there are more photons and neutrinos 
generated by the particle decay than electrons and protons. In the regime 0.2 < x < 0.4, photons, 
neutrinos and protons are generated with roughly equal abundances. However, for x <^ 0.5, photons 
and neutrinos again outnumber protons, in particular protons cut off at x ~ 0.75 whereas neutrinos 
and photons are generated in the cascades with energies up to x ~ 0.95. These differences will lead 
to different shapes of the expected fluxes h{E) as can be seen from eq.(12). 

We now compare our proton fragmentation function with the commonly used Hill approxima- 
tion in Figure |j. Although his form provides a good fit for a low decaying particle mass, viz. 
mx = 10 3 GeV, it no longer does so for a high mass, viz. mx = 10 11 GeV. This is understandable 
given that the numerical co-efficients in eq. ( 12 ) were chosen to match relatively low energy collider 
data. However the functional form itself is well motivated and using our HERWIG runs we can 
determine new numerical co-efficients appropriate to heavier mass particles. Another advantage of 
the present approach is that the spectrum of neutrinos and photons is determined separately from 
that of the protons and not simply assumed to be proportional as in previous work OH]. 

To study the highest energy cosmic ray events we need to consider particle masses beyond 
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10 11 GeV but this is difficult to do directly with HERWIG for technical reasons as mentioned 
earlier. We therefore resort to an extrapolation procedure as follows. For the range x G [0, 0.2] 
the fragmentation functions are smooth and evolve monotonically with mx so the fragmentation 
functions for a 10 13 GeV particle is obtained from simple linear extrapolation of the lower energy 
fragmentation functions in each individual energy bin. For x £ [0.2, 0.6] we first fit the calculated 
fragmentation functions to the form 



for protons, and the form 



diVfiti = exp[c 2 yin(l/x7](l-x) 2 
dx Cl xy/]n(l/x) ' 1 ' 



diVfit2 = d exp[rf 2 ln(l/x)](l - xf 
dx 1 xln(l/x) 



which proves more suitable for photons, neutrinos and electrons. The numerical co-efficients cj, c 2 
and d\ , d<i are determined for particle masses less than 10 11 GeV by minimizing \ 2 m the fit to 
the actual HERWIG runs. In Figures | and | we show these fits for x < 0.6 to the proton and 
neutrino fragmentation functions corresponding to masses of 10 5 and 10 9 GeV. Then we determine 
the appropriate co-efficients for heavier masses by extrapolation. An example, for mx = 10 13 GeV, 
is shown in the figures. For x ^ 0.6, statistical fluctuations become too severe so we extrapolate 
the fitting functions between the value at x = 0.55 and a cutoff which is taken to be x = 0.75 for 
protons and x = 0.95 for neutrinos, based on the observed behaviour for masses upto 10 11 GeV 
shown in Figures || and ||. 

Finally, we mention the continuation of the proton fragmentation function for very low x, viz. 
x < 0.015, which is relevant at high masses e.g. mx = 10 13 GeV. Since it proved impractical to 
have additional binning intervals at very small x, we employ the fragmentation function (|l4|) in 
this regime, normalized to our computations at x = 0.015. 



B. Comparison with observations 

With the fragmentation functions obtained above, we can now calculate the expected fluxes 
of protons and neutrinos from decays of particles such as cryptons in the halo. We normalize the 
calculated proton flux fll2] ) to the observed cosmic ray flux at 10 19 eV : 

log 10 [ip(10 19 eV) /m~ 2 sec^sr" 1 eV 2 ] = 24.32. (18) 

Note that the corresponding neutrino flux I V {E) is then a prediction as the fragmentation function 
for neutrinos is computed independently. 

The expected proton fluxes are shown in Figure @. We see that a crypton with mx = 10 11 GeV 
fits the flat power law well but cannot explain the events beyond 4 x 10 19 eV. Although this is easily 
achieved for mx = 10 13 GeV, the decays of such a massive particle would overproduce protons for 
E > 3 x 10 19 eV. Thus a crypton with mass mx ~ 10 12 GeV provides the best compromise although 
it too predicts a spectrum somewhat flatter than the one indicated observationally. (The reader is 
reminded that all differential fluxes have been multiplied by E 3 in eq.fll^).) 

An interesting signature for forthcoming experiments is the predicted ratio of the proton 



to neutrino flux [70|. In Figure E we compare the expected flux of protons and neutrinos for 



10 



mx = 10 GeV. (We also show the photon flux to illustrate the difference from the predic- 



tion in ref. [37] but emphasize that this will be degraded through interactions with photon back- 
grounds during travel to Earth.) As can be seen, the neutrino flux exceeds the proton flux for 
10 19 eV < E < 10 20 eV and also for E > 3 x 10 20 eV, as may have been anticipated from the 
comparison of their respective fragmentation functions. Thus the ratio I p /I u has a characteristic 
peak at about 2 x 10 20 eV as shown in Figure ||. This could be an useful diagonistic of the decaying 
particle hypothesis for future experiments such as the Pierre Auger Project. Note that taking the 
extragalactic contribution into account would boost the neutrino flux by a factor of ~ 2 over that 
shown in the figures. 

The abundance and lifetime of decaying particles such as cryptons are related through the 



spectrum normalization ( |18[) as: 

log 10 (/cos^x/i 2 ) = k + log 10 (^j (19) 

where k = —15.78, —16.13, —15.59 for crypton masses mx = 10 11 , 10 12 , 10 13 GeV respectively. For 
a given crypton mass, a higher lifetime must be compensated for by a higher relic abundance, as 



illustrated in Figure 10. So for example, if / CO s^x h 2 ~ 1, cryptons with a mass of mx = 10 12 GeV 
are required to have a lifetime of tx ~ 10 16 yr if they are to explain the EHECR flux. If the 
enhancement in the halo is / cos ~ 3 x 10 4 as expected for cold dark matter, then the lifetime may 
be increased to ~ 4 x 10 20 yr if fix h ~ 1; alternatively, for the same lifetime one could tolerate a 
lower relic abundance Qx ^ 2 ~3x 10~ 5 . 

With regard to the fluxes of electrons and photons, both species would generate electromagnetic 
cascades on the prevalent radiation backgrounds through pair production and inverse Compton- 
scattering. A thorough analysis of such propagation effects and the resulting modifications of the 



injected photon and electron spectra has been performed [7_1|. It was found that the relic decaying 
particles with mx ^ 10 14 GeV would contribute excessively to the diffuse 7-ray background and 
are therefore ruled out. Hence, the mass range we favour, viz. 10 11 < mx/ GeV < 10 13 , does not 
lead to any conflict with observations. This conclusion is strengthened by the fact that according 
to our calculations the previous estimate |37j of the 7-ray flux from decaying particles was too 
high. Although the positrons released in the decays may be accumulated in the galactic halo, the 
astrophysical uncertainty in the containment time does not allow a restrictive constraint to be 
derived from limits on the positron flux in cosmic rays [f72^ [. 

With regard to the neutrino background, the predicted flux at high energies is well below the 
upper limits derived from consideration of horizontal air showers [ 73 1 , again because the decaying 



crypton mass is restricted to be less than about 10 13 GeV. It is also interesting to consider the 
flux at lower energies of O(10 3 ) GeV where experiments such the forthcoming ANTARES detector 



[74] will be most sensitive. As seen in Figure || the predicted neutrino flux dominates over the 
proton flux at low energies, thus the bulk of the energy released by the decaying cryptons ends 
up as neutrinos. ^ Therefore we expect the neutrino flux at TeV energies to be at least ~ 10 8 



4 This expectation motivated the study undertaken earlier in which the abundance/lifetime of 
massive metastable relic particles was constrained using experimental limits on the high energy 
neutrino flux set by underground nucleon decay detectors and the Fly's Eye experiment |5C|. 
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times larger than the EHECR flux at 10 eV. Moreover the neutrinos should be well correlated 
in both time and arrival direction with the cosmic rays since the path length in the galactic halo 
is < 100 kpc. This is in contrast to the case of other suggested cosmologically distant sources such 
as gamma-ray bursts where the relative time delay can be upto ~ 10 3 yr |ff5| . 



IV. CONCLUSIONS 

We have investigated the hypothesis that the highest energy cosmic rays, in particular those 
observed beyond the GZK cutoff, arise from the decay of massive metastable relic particles which 
constitute a fraction of the dark matter in the galactic halo. To simplify computations (using 
the HERWIG Monte Carlo event generator) we have considered only decays into qq pairs with 
unit branching ratio. Comparison with experimental data indicates that a decaying particle mass 
of 0(1O 12 ) GeV is required to fit the spectral shape while the absolute flux requires a lifetime of 
0(lO 2O )yr if such particles contribute the critical density. The predicted decay spectra may be 
somewhat altered if 3-body decays and other final states (e.g. supersymmetric particles (7(||) are 
considered. However our conclusions regarding the preferred mass and relic abundance/lifetime of 
the decaying particle are unlikely to be affected. In particular it would appear that the approxima- 



tions used to calculate the particle spectra in previous studies of decaying topological defects [34] 



and hypothetical massive particles [37| were not sufficiently accurate. Our work indicates that the 



topological defect model is disfavoured unless the mass of the decaying gauge bosons is less than 
about 10 13 GeV, which is well below the unification scale of ~ 10 16 GeV. (A similar conclusion is 
arrived at by independent arguments in refs. |77 ,ff8|.) By contrast, cryptons from the hidden sector 



of supersymmetry breaking have a mass of the required order, as well as a decay lifetime which is 
naturally suppressed. However their relic abundance is difficult to estimate reliably, although we 
have argued that it may be cosmologically interesting. 

The primary intention of this work is to attempt to quantify the decaying particle hypothesis 
in a manner which is of interest to experimentalists. We have therefore computed the expected 
neutrino to proton ratio as a function of energy since this is an important test of competing 



hypotheses for forthcoming experiments, in particular the Pierre Auger project [69|. Of course 
our cleanest prediction is that the cosmic ray spectrum should cut off just below the mass of the 
decaying crypton, at ~ 3 x 10 20 eV. Moreover, with sufficient event statistics it should be possible 
to identify the small anisotropy which should result from the distribution of the decaying particles 
in the Galactic halo J7^]. Thus although the hypothesis investigated here is very speculative, it 
is nevertheless testable. Perhaps Nature has indeed been kind to us and provided a spectacular 
cosmic signature of physics well beyond the Standard Model. 
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FIG. 1. The high energy cosmic ray spectrum beyond the 'ankle'. (Note that the differential 
flux has been multiplied by E s .) The data shown are from AGASA, stereo Fly's Eye, Haverah 
Park and Yakutsk and are AGASA-normalized ||. The highest energy monocular Fly's Eye event 
at 3 x 10 19 eV is also shown. A fit to the spectrum from the superposition of a steeply falling and 
a flatter power law (dashed lines) is indicated Q. 
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FIG. 3. Fragmentation functions of photons, neutrinos, and electrons compared to that of 
protons for a decaying particle mass of mx = W n GeV. 
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FIG. 4. Comparison of the computed proton fragmentation function with the leading-log ap- 
proximation of eq. (p^ ), normalized to the HERWIG computation at x = 0.042. The upper solid 
line refers to a decaying particle mass of mx = 10 3 GeV and the lower one to mx = 10 11 GeV. 
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FIG. 5. Extrapolation of the computed proton fragmentation function to higher decaying parti- 
cle masses. The upper two solid lines are HERWIG results for decaying particles of mass rax = 10 5 
and 10 9 GeV while the dashed lines are the best fitting functions according to eq. (|l6|) . (The cor- 
responding x 2 -values are also indicated.) The lower solid line is the fragmentation function for 
mx = 10 13 GeV obtained from extrapolating the fitting parameters. 
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FIG. 7. Expected proton flux for decaying particle masses mx = 10 11 , 10 12 and 10 13 GeV com- 
pared with observations. The theoretical spectra are normalized at 10 19 eV to the flat component 
(dashed line) suggested by the Fly's Eye data. 
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FIG. 8. Expected fluxes of protons, neutrinos and photons from decays of a decaying particle 
with mass mx = 10 12 GeV (normalized as in Figure 0) compared with the observations. Note that 
the photon flux will be degraded through interactions with the CMB during travel to Earth and is 
shown for illustrative purposes only. 
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FIG. 9. The ratio of the proton to neutrino flux for the decaying particle mass mx = 10 12 GeV. 
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FIG. 10. The relic decaying particle abundance versus lifetime for various masses, as required 
by the flux normalization to the observations. 



25 



